Action Potential Onset Dynamics and the Response Speed of Neuronal Populations 



B. Naundorf, T. Geisel and F. Wolf 
Max- Planck- Institut fur Strdmungsforschung and Fakultat fiir Physik, 
Universitat Gottingen, 37073 Gottingen, German^ 
(Dated: 9th February 2008) 

The result of computational operations performed at the single cell level are coded into sequences 
of action potentials (APs). In the cerebral cortex, due to its columnar organization, large number of 
neurons are involved in any individual processing task. It is therefore important to understand how 
the properties of coding at the level of neuronal populations are determined by the dynamics of single 
neuron AP generation. Here we analyze how the AP generating mechanism determines the speed 
with which an ensemble of neurons can represent transient stochastic input signals. We analyze a 
generalization of the 6-neuron, the normal form of the dynamics of Type-I excitable membranes. 
Using a novel sparse matrix representation of the Fokker-Planck equation, which describes the 
ensemble dynamics, we calculate the transmission functions for small modulations of the mean 
current and noise noise amplitude. In the high-frequency limit the transmission function decays 
as OJ where 7 surprisingly depends on the phase 9s at which APs are emitted. If at 9s the 
dynamics is insensitive to external inputs, the transmission function decays as (i) lu~^ for the case 
of a modulation of a white noise input and as (ii) uj~'^ for a modulation of the mean input current 
in the presence of a correlated and uncorrelated noise as well as (iii) in the case of a modulated 
amplitude of a correlated noise input. If the insensitivity condition is lifted, the transmission function 
always decays as <^~^, as in conductance based neuron models. In a physiologically plausible regime 
up to IkHz the typical response speed is, however, independent of the high-frequency limit and is 
set by the rapidness of the AP onset, as revealed by the full transmission function. In this regime 
modulations of the noise amplitude can be transmitted faithfully up to much higher frequencies 
than modulations in the mean input current. We finally show that the linear response approach 
used is valid for a large regime of stimulus amplitudes. 



I. INTRODUCTION 

Neurons are the basic building blocks of neural net- 
works and thus constitute the computational units of the 
brain. They dynamically transform synaptic inputs into 
output action potential (AP) sequences. To conceive the 
complex computational capabilities of the brain, it is cru- 
cial to understand this transformation and to identify 
simple neuron models which accurately reproduce the dy- 
namical features of cortical neurons. 

Here we study this mapping in a reduced neuron 
model. This model is obtained by a generalization of 
the 6'-neuron fill. IT3|. which is a canonical phase oscil- 
lator model of excitable neuronal membranes exhibiting 
Type-I excitability. Phase oscillator models have a long 
history in physics and biology 0, 0, [i^, and re- 
cently they were introduced in theoretical neuroscience 
[Tll |. In contrast to integrate- and-fire models, which are 
phenomenological models of cortical neurons, they can 
be derived from the limit cycle dynamics of conductance 
based neuron models, reducing the complex dynamics 
which usually incorporates many degrees of freedom to 
a single phase variable. This reduction is an important 
prerequisite for analytical studies of either the dynamics 
of single neurons or of neural networks. 

Cortical neurons in vivo are subject to an immense 
synaptic bombardment, resulting in large fluctuations of 
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their membrane potential (MP) 0, and irregular 

action potential flring j^^. Because the exact compu- 
tational role of these fluctuations is largely unknown, it 
was suggested to treat them as a random process, divid- 
ing the synaptic input into a mean input current and a 
random fluctuating contribution with a given cor- 
relation time Tc. The fluctuations can serve as a po- 
tentially independent information channel because when 
the afferent activity of a neuron changes, not only the 
mean input is affected, but also the amplitude of the 
fluctuations jlElliH. 

The stationary response properties of the classical 6- 
neuron subject to fluctuating input currents were calcu- 
lated in (l§ . l2a | for a temporally uncorrelated input and 
in 01 for a temporally correlated input current. Both 
studies showed that the 6'-neuron can reproduce the sta- 
tionary response properties exhibited by many cortical 
neurons, i.e. a square-root dependence of the flring rate 
on the input current close to threshold for small noise 
amplitudes and irregular firing in the noise driven 
regime. Despite its success to reproduce the stationary 
firing behavior of cortical neurons, the 6'-neuron lacks a 
crucial dynamical feature: The fast action potential up- 
stroke exhibited by conductance based neuron models. 
Here we study a generalization of the classical 6'-neuron 
with an adjustable action potential onset speed, intro- 
ducing a phenomenological term which mimics the fast 
activation of sodium channels. 

We derive the time dependent response in the pres- 
ence of temporally correlated noise to a modulation in 
the mean input current and a modulation in the noise 
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amplitude. For both modulation paradigms we calcu- 
late the high frequency limit. In this limit, the response 
amplitude decays as where the integer exponent 7 

is completely independent of the action potential onset 
dynamics and surprisingly only depends on the oscilla- 
tor phase 9s^ at which an action potential is emitted: If 
at 9s the dynamics is insensitive to external inputs, the 
transmission function decays as (i) ijO~^ for the case of 
a modulation of an uncorrelated noise amplitude and as 
(ii) [jj^^ for a modulation of the mean input current in 
the presence of a correlated and uncorrelated noise as 
well as (iii) in the case of a modulated amplitude of a 
correlated noise input. If the insensitivity condition is 
lifted, the transmission function always decays as 
as in conductance based neuron models. 

The full transmission function is then calculated via 
the eigenvalues and eigenfunctions of the Fokker-Planck 
operator, which describes the dynamics of the ensem- 
ble averaged probability density function. The eigenval- 
ues and eigenfunctions are computed using a high per- 
formance iterative scheme, the Arnoldi method |2l. 
from a novel sparse matrix representation of the Fokker- 
Planck operator. This method allows for a fast compu- 
tation and high numerical precision, hard to achieve by 
direct numerical simulations. 

We then demonstrate that the response amplitudes for 
the classical 6'-neuron typically exhibit a cut-off behav- 
ior, where the cut-off frequency, which is closely linked 
to the spectral properties of the Fokker-Planck operator, 
is approximately given by the neurons stationary firing 
rate. Stimulations at frequencies larger than the cut- 
off frequency are strongly damped. For an increasing 
action potential onset speed at a fixed stationary rate, 
stimuli with much larger frequencies can be transmitted 
almost unattenuated. We show that the response ampli- 
tude for the case of a noise modulation typically decays 
much slower than in the case of a mean input current 
modulation. 

The impact of noise on the dynamic response prop- 
erties was previously almost exclusively studied in 
integrate-and-fire models [IS ■ The first studies were 
pioneered by Knight who considered a simple in- 
tegrator model, in which the firing threshold is drawn 
randomly, every time an action potential occurs. These 
results were then extended to models, where the reset 
voltage was also drawn randomly, and to models in which 
the the input changed either very slowly, or to spike re- 
sponse models, where the input is assumed to change very 
fast Recently, the impact of current noise on the 

dynamical response of the leaky integrate-and-fire model 
was investigated 0, HjQ^SES- these studies it was 
shown that integrate-and-fire models driven by a synap- 
tic fiuctuating input exhibit a linear response amplitude 
which does not decay to zero in the high frequency limit. 
This lead some to the conclusion that cortical neurons 
can transmit information instantaneously 0, ^23] . Only 
recentlv j th is interpretation was questioned by two stud- 
ies [1^ Efl| which demonstrated that the unattenuated 
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Figure 1: Phase plane of a type-I single compartment conduc- 
tance based model (Morris-Lecar model 12.8] ) in the excitable 
regime (filled dot: stable fixed point, open dots: unstable 
fixed points). Gray lines are the nullclines, denoted by it) = 
and y = 0. Black lines are stable and unstable manifolds 
of the saddle and the node. The excitable dynamics can be 
reduced to a phase oscillator with one degree of freedom pa- 
rameterized by the angle 

transmission of high frequency signals in integrate-and- 
fire models are more a consequence of the oversimplifica- 
tion of the model rather than property of real neurons. 

II. MATERIAL AND METHODS 
A. Model 

The model is based on the normal form of the dynam- 
ics of type-I membranes at the bifurcation to repetitive 
firing. Conductance based neuron models which exhibit 
Type-I excitability typically undergo a saddle-node bi- 
furcation of codimension one, when brought to repetitive 
firing. A center-manifold reduction at the bifurcation 
point leads to the following normal form (3^ : 

CT> = -^(^-l^*)' + (/W-^c), (1) 

which is a dynamical equation for the MP V of the neu- 
ron. The input current relative to the rheobase Ic of the 
neuron is denoted by lit). The constants A and V* can 
be deduced from a given multidimensional conductance 
based model. It is convenient to introduce dimensionless 
quantities Vand /: 

V = {y-v*)iv^ (2) 

m ^ {m - / (gVo) (3) 
and the effective time constant: 

r = C/g (4) 
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The rescaled dynamics is then given by: 

tV ^V'^ + I{t) (5) 

For > 0, the MP has a finite "blow-up" time, mean- 
ing that it needs a finite time to get from — oo to +oo, 
where both ends of the real axis are identified, turning the 
model into a phase oscillator. The normal form Eg . O 
is equivalent to a phase oscillator, the 6'-neuron ll^ . 
Its equation of motion, 

r^= (l-cos6l)+/(i)(l + cos6') (6) 

is found by substituting V = tan(6'/2) with the angle 
variable 9 in the interval (— 7r,7r]. 

In the model, a spike is emitted each time 9 reaches 
the value 9s. By choosing 9s = n, the original 0-neuron is 
obtained. Figurenjillustrates schematically the reduction 
of a conductance based neuron model to a phase oscillator 
model. 

Although the 6'-neuron is the normal form of the dy- 
namics at the bifurcation, it lacks the rapid AP onset 
exhibited by conductance based neuron models and real 
neurons. To account for this dynamical feature we gen- 
eralized the model to refiect the rapid depolarization of 
the membrane resulting from the fast kinetics of sodium 
conductances in the following way: 

tV = V^ +Iit) + a{l + taiih{(3V)), (7) 

where we introduced two additional parameters a and 
p. The sigmoidal term phenomenologically models the 
part of the sodium channel activation curve, which is 
not included in the y^-term of the normal form. The 
parameter a controls the sodium peak conductance and 
the parameter (3 the width of this activation curve. Both 
parameters control the rapidness of the AP onset. As for 
the normal form an equivalent phase oscillator equation 
can be found by substituting V = tan {9/2): 

t9 = (1 - cos6') + (1 + cos6i) • (8) 
-I- {I{t) + a (1 + tanh(/3 tan (61/2))} 

B. Fluctuating input currents 

In vivo., neurons are subject to an ongoing synaptic 
bombardment, resulting in a fiuctuating MP. To model 
this situation, we assume a temporally fiuctuating input 
current, 

I{t) = h + az{t), (9) 

composed of a mean /q and a stationary fiuctuating part 
az{t), where z{t) is an Ornstein-Uhlenbeck process with 
a given correlation function {z{t)z(t')) = exp(— i/rc). 
Thus z{t) obeys the Langevin equation 

TAz{t)^-Z+V^Tl{t) (10) 



with {'q{t)) = and {'q(t)r](t')) = S(t - t'). Eq. © and 
Eq. {HI describe a realization of the dynamics of a sin- 
gle neuron. Since the input is fiuctuating and we are 
interested in coding at population level it is natural to 
consider an ensemble of such units, described by the time 
dependent probability density function P{9,z,t). Its dy- 
namics is determined by the Fokker-Planck equation |32 |: 

dtP{9,z,t) = LP{9,z,t), (11) 

with, 

L = -T-^de{{l-coa9) 

+ {lo + <7z + a{l + tanh(/3tan(6'/2)))) 

■ {l + coa9)} + T~^d,(^z+^d,y (12) 

The boundary conditions for P{9, z, t) are periodic in the 
9- and natural in the z-direction. 



C. Time dependent firing rate 

The ensemble averaged firing rate is given by the prob- 
ability current across the line 9 — 9s with positive veloc- 
ity. At = TT the dynamics is independent of the input 
current I[t) and the rate is equal to the probability cur- 
rent through the entire line 9 — it: 

/oo 
dzP{n,z,t) (13) 
-oo 

Although quite convenient for analytical considerations, 
the definition of this spike-phase is, however, rather arbi- 
trary. In the normal form, the point 9s = tt corresponds 
to the point V = oo, where the model refiects least the 
dynamics at the bifurcation. To assess if this particu- 
lar choice has any infiuence on the dynamical response 
properties of the model, we also calculate the firing rate 
at 9s = TT — S. The probability current through this fine 
is given by: 

/oo 
P{9s,z,t){{l-cos9s) (14) 
-OO 

+ (/o + (Tz + a(l + tanh(/3tan(6's/2)))) dz 

The rate is, however, not exactly given by the fiux 
Jg. There is a contribution from trajectories, which are 
driven back below the threshold due to the external fiuc- 
tuations. For a correlated input current, however, the 
introduced error is exponentially small. This can be seen 
in Eq. JHl- For small values of S, the probability distri- 
bution P{9, 9) around tt — 5 is a Gaussian with a mean 
value 2 — 5^ and a width oc 5^ . The negative part of this 
Gaussian is proportional to: 
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For all practical purposes {S < 0.5 and cr < 1), this inte- 
gral is smaller than 10~^°. We will see, however, that the 
definition 9s — tt in the classical ^-neuron qualitatively 
changes the dynamic response of the model in the high 
frequency limit. 



D. Parameter choice 

Before discussing the stationary and dynamical proper- 
ties of the generaHzed 6'-neuron we would like to define a 
biologically plausible parameter regime. The parameters 
which we need to fix are the time constant r, the mean 
input current /q, the strength of the fiuctuating input a 
and the synaptic input correlation time Tc. An estimate 
of the correlation time of the MP is given by approximat- 
ing the dynamics for Iq < near the stable fixed point 
by an Ornstein-Uhlenbeck process. Straightforward lin- 
earization around the stable fixed point at — VTo then 
yields: 



'T'rclax 



(16) 



In the subthreshold noise-driven regime, which we will 
discuss in the following, we choose Iq = —0.1. The time 
constant r is then adapted via Eq. ifTfijl . to achieve a re- 
laxation time of approximately 5ms, which leads to values 
for T of approximately 3ms. 

The parameters a and /? parameterize the sodium ac- 
tivation curve, which, in conductance based models, de- 
termines the speed at the action potential onset. For the 
following numerical treatment we keep /?, which mediates 
the width of the activation curve and is an intrinsic phys- 
iological parameter, fixed to a value of 20. The parame- 
ter a, which represents the sodium peak conductance, is 
changed in the range from to 1. 

Figure |2l shows three sample realizations of Eqs. (jHl 
EJ for different values of the parameter a. If the input 
current is positive for a sufficient amount of time, action 
potentials are initiated. With increasing values of a the 
sharpness at the onset increases, while the subthreshold 
fiuctuations are not affected. 



E. Dynamic Response Theory 

For time-dependent input currents el{t), the Fokker- 
Planck operator L{9,z,t) can always be split into two 
parts: 



L{0,z,t) = Lo{9,z)+eLi{e,z,t), 



(17) 



where Lo{0,z) is the time- independent part and 
Li{0,z,t) contains all time-dependencies of the exter- 
nal input. In the following we require that the time- 
dependent inputs are small in magnitude, i.e. e ^ 1. 
We then expand the general time-dependent solution in 




Time (ms) 

Figure 2: Increasing a leads to a sharper action potential 
onset, (a) Sample MP trajectories for a = 0, a = 0.1 and 
a — 1. The inset shows the deterministic part of Eq. 
(b) Fluctuating input current I{t). Parameters: Tc = 1.5ms, 
a = 0.3, Jo = -0.1 and /3 = 20. Right before AP onsets the 
trajectories are virtually identical. 



powers of e : 

PTDiO,z,t) = Pa{e,z)+ePi0,z,t)+Oie'') (18) 

Inserting this solution into the Fokker-Planck equation 
and keeping only terms up to linear order in e leads to 
a dynamical equation for the time dependent part of the 
density P{0, z, t): 

dtPi0,z,t) = Lo{9,z)P{e,z,t) + Li{9,z,t)Pa{e,z) (19) 

Formally the solution of this equation is given by: 



P{0,z,t) 



Li{0,z,t)Poi0,z)dt' (20) 



In the following we will consider stimuli of the type: 

Li{e,z,t) = e"^'Li{0,z) (21) 
Eq. lf2(Hl can then be readily solved, yielding: 



Pi0,z,t)^j2 



Ck 



iuj — A; 



-Pk{0,z)e'' 



(22) 



The Ck are the expansion coefficients of Li{6, z)Po{0, z) 
into the eigenfunctions Pk{0,z) of Lq{9,z). The time- 
dependent firing rate is given by Eq. lfT3ll : 



u{t) 



/OO 
dz ((l-cos6's) 
-OO 
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Figure 3: Sketch of the population response paradigm. An ensemble of neurons receives a modulated noisy input current or 
a current, where the noise amplitude is modulated. The noise realization which each neuron receives is different, leading to 
different MP traces and AP sequences. The output quantity is the population averaged firing rate in the interval [t,t + dt), 
u{t). 



+ {Iq + az + a{l + tanh(/3 tan(6's/2))))) 

=: z/o + £i^i(c^)e'('^*+'^('^» (23) 

In the following we will consider two types of external 
stimulations: 

1. Modulations in the mean input current: 

Jo ^ /o + ee^"' 

2. Modulations in the noise amplitude: 



F. High frequency limit 

In this section we sketch how to analytically calculate 
the asymptotic decay of h'i{u!) in the limit uj oo. In- 
serting Eq. (|22l into Eq. leads to: 

(^iuj - Lo) P{0, z, t)e'"^' = LiPo{9, z) (24) 

If the right hand side vanishes at 6* = ^?s, P{9,z,t) has 
to decay at least as uj~'^. Differentiation of Eq. with 
respect to t and subsequent reinsertion leads to: 

+ Lo) P(0, z, t) = -LoLiPo{9, z) (25) 

If now the right hand side vanishes at ^? = 6*^, Eq. Ijl9|l 
has to be differentiated again, until, after reinsertion, the 
right hand side is different from zero. 

G. Matrix Method 

As demonstrated the dynamical response properties of 
the generalized 6'-neuron to small time-dependent inputs 



are completely determined by the spectrum and eigen- 
functions of the Fokker-Planck operator L. To compute 
the dynamical response properties in the presence of a 
temporally correlated noise current for arbitrary stimu- 
lation frequencies we expand L into a complete orthonor- 
mal basis leading to a sparse matrix representation for 
which we compute the eigenvalues and eigenfunctions nu- 
merically. This approach has the advantage that the re- 
sponse properties can be computed with very high accu- 
racy. The two subtleties we will have to deal with are 
that (1) the resulting matrix is very large in the parame- 
ter regime we are interested in (up to 10^ x 10^) and (2) 
the operator L is not Hermitian and thus standard di- 
agonalization procedures such as the Lanczos algorithm 
can not be applied. We solved both problems by using a 
basis-set, which results in a very sparse matrix represen- 
tation, and by using a high performance iterative scheme, 
the Arnoldi method to compute the eigenfunctions 
and the spectrum of this matrix to a high numerical ac- 
curacy. 



H. Eigenvalues and eigenfunctions for a correlated 
noise input 

1. Matrix equation 

We first replace the probability density P{9, z, t) in an 
eigenmode Ansatz with e^''*^Pk{9, z). Inserting this into 
Eq. the exponential cancels: 

XkPk{9,z) = LPk{9,z) (26) 

Due to the imposed boundary conditions, the set {Afc}, 
i.e. the spectrum of L{9, z) is discrete. There is, however, 
a macroscopic drift in the system, meaning that detailed 
balance is not fulfilled and thus L is not Hermitian [l5| . 
This means that the resulting spectrum {A^} and the 
corresponding eigenfunctions Pk{9,z) are complex. By 
complex conjugation of Eq. Ij26|l it is easy to show that 
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to every eigenvalue Xk with the corresponding eigenfunc- 
tion Pk{9,z), an eigenvalue with the eigenfunction 
P^{9,z) exists. This guarantees that a real solution can 
always be constructed. The solution with Aq = corre- 
sponds to the stationary density and the time dependent 
solution can always be given in terms of eigenfunctions 
and eigenvalues 



(27) 



with Pinitiai(6', 2;) = J2k°'kPk{d,z). Although the eigen- 
functions of L form a basis, it is important to note that 
they are not orthogonal. An important property is that 
the mean value of all eigenfunctions except Po{6,z) is 
zero: 



dd / dzPk{e,z) = 

-TT J — 00 



(28) 



To actually compute the spectrum and eigenfunctions we 
expand P{9, z) into a set of complete orthonormal func- 
tions: 



00 

P{6,z) = ^ an^yn'>i^n,m{d,z) 
m=0 



(29) 



with 



in9 TT 



-1/2 



(30) 



This expansion obeys the imposed boundary conditions. 
In the 6'-directions it consists of plane waves, while in 
the z-direction harmonic-oscillator functions are used 
with the Hermite polynomials Hmiz) f^. We now in- 
sert Eq. (|2ni in Eq. (tTT|l . Multiplying from left with 
V'n' m' ^nd integrating over the whole domain leads 
to a matrix eigenvalue equation for the (a„^m): 

Aa„,,„ = {~iT'^{l + Io)n~T~^m)an,m 

+ (2r)"-^i(l - Io)n (a„„i,m -|- an+i,m) 
ina , 

[m + l)a„^,„+i + ma-n^m-i 



1 



+ -{m + 1) (a„_i,„i+i -I- an+i,m+i) 
+ -{m + 1) (a„_i,m_i -I- an+l^m-l) 



+ {V2tTc) ^ V(m + l)(TO + 2)a 



n,m+2 



-ZT a Qr, 



1 ^ 



k^l 



{iCk — Sk) au+k,T. 



/ ^ -^n,m;n' ,m' ^n' ,m' • 
7i' ,7n' 



(31) 



The coefficients C^and Sk denote the Fourier compo- 
nents of (1 -I- COS0) tanh (/3tan {9/2)) of the expansion in 
cos {k9) and sm{k9) up to order K. To solve this eigen- 
value problem numerically we have to restrict the indices 
n and m to 



n£ {~N ...N}, me{O...M} 



(32) 



Since the stationary density is very peaked for realistic 
firing rates, we need many plane wave basis functions, 
i.e. up to iV « 10^. With M = 50 the matrix that we 
have to diagonalize will be of size 10^ x 10^. To only rep- 
resent this matrix in full form would require 3.8- lO'^GB of 
storage capacity. We note, however, that the matrix L in 
Eq. Ij31|l is very sparse, for a = it connects an element 
an,m even only to the elements a„±i,m±iand an,m+2- For 
a > the number of nonzero entries in L solely depends 
on the number of Fourier components K of the AP on- 
set term of the generalized model. In general, however, 
the number of elements in the matrix L is only of order 
NxM, i.e. very sparse compared to its full size x M^. 
This makes it possible to use a high performance itera- 
tive algorithm, the Arnoldi-method 0, to solve this 
eigenvalue problem numerically. The time-dependent fir- 
ing rate ^{t) is calculated using Eq. Ij23|l . 



III. RESULTS 
A. High frequency limit 

1. Dynamics insensitive at action potential (9a = n) 

For both types of input modulations the modulus of the 
right hand side of Eq. Ij24|l vanishes at 6* = tt. Therefore 
the P{9,z,t) has to be at least of order uj^'^, such that 
the left hand side vanishes for uj ^ 00. Differentiation of 
Eq. Ijl9|l and subsequent reinsertion leads to: 



Lo]Pi9,z,t) = -LoLiP„i 



(33) 



The right hand side does not vanish at 9 — tt in the 
case of a mean current modulation and in the case of a 
modulation in the noise amplitude. Since both sides have 
to be real valued, the modulus of P{9, z) has to be cx uj^'^ 
and the phase </3(w) goes to —it. 

In the limit Tc 0, i.e. an uncorrelated input current, 
the same argument holds in the case of a mean current 
modulation. For a modulation in the noise amplitude, 
the right hand side of Eq. Ij33|l is zero, resulting even in 
a decay and a phase lag of 37r/2. 



2. Generic case (9s 7^ tt) 

For 9s ^ TT - S, S > the right hand side of Eq. lf2ijl 
does not vanish. This means that for large frequencies 
the rate modulation 1^1(00) decays as w"^ and the relative 
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Table I: High frequency behavior of the generalized 0-neuron, 
the leaky integrate-and-fire model 0, and conductance 
based models. The response of a conductance based model 
for Tc > and a mean current modulation was studied in [T3 |. 
The asymptotic response of the conductance based model in 
the other cases follows from the same argument as for the 
asymptotic response of the S-neuron and was confirmed by 
direct numerical simulations (data not shown) 



phase shift ipito) is — 7r/2, which is the same asymptotic 
decay as in conductance based neuron models. Table 
in summarizes the high frequency behavior of the gener- 
alized 6'-neuron and compares it to the high-frequency 
limit of conductance based model neurons and the clas- 
sical leaky integrate-and-fire model. 

We would like to point out, that the cj"^ and iu~^ 
decay of the classical 6'-neuron is only due to (i) the in- 
sensitivity of the dynamics to inputs at 9 = n and the 
symmetric up- and downstroke of the action potential 
around 6s = tt. Here, both conditions are lifted by defin- 
ing the spike phase at a different value than tt. Another 
way to induce a w~^-decay would be to change the right 
hand side of Eq. l|Sl, such that LiPq does not vanish at 
9 = n, e.g. by introducing high order terms in cos6'. 
This would however require a structural change of the 
oscillator dynamics. A second important point to note 
is the independence of the high-frequency limit from the 
dynamics at the action potential onset. 



B. Linear response transmission Functions 

Using the matrix method described above, we com- 
puted the linear responses to modulations in the mean 
input current and to modulations in the noise amplitude. 
Figure 31 summarizes the response amplitude curves for 
the 0-neuron model, the generalized 6'-neuron model and 
compares them to direct numerical simulations of the re- 
sponse of the leaky integrate-and-fire (LIF) model. 

The 6'-neuron exhibits a cut-off behavior in its response 
amplitude to both types of input modulations. Frequen- 
cies up to the stationary firing rate can be transmit- 
ted unattenuated larger frequencies are strongly damped. 
For an increasing onset speed and fixed stationary rate 
the resonance maximum shifts only to slightly larger fre- 



quencies, a dramatic change, however, occurs at interme- 
diate frequencies up to IkHz. In this regime the response 
amplitude is substantially lifted to much larger transmis- 
sion amplitudes. This effect is much more pronounced for 
the case of a modulation in the noise amplitude than for 
modulations in the mean input current, leading to an 
undamped response for frequencies up to 200Hz. The 
LIF model, on the other hand, shows a completely arti- 
ficial response behavior. The transmission function, for 
both types of modulations does not decay at all, even 
for frequencies up to IkHz. For modulations in the noise 
amplitude, the transmission function can even grow for 
increasing stimulation frequencies. 



C. Nonlinear response for large stimulation 
amplitudes 

So far we have only considered the Hnear response 
transmission function, which is strictly speaking only 
valid in the limit in which the stimulation amplitude goes 
to zero. Here we show, however, that the linear response 
covers a large range of input amplitudes. In principle, 
we could use the same matrix method employed for the 
linear response theory, taking into account higher order 
Floquet modes Here we explore this regime, how- 

ever, by direct numerical simulation of Eq. Figure 
El shows the amplitude of the first four Fourier modes of 
the rate response as a function of the overall amplitude of 
the rate modulation. For both types of modulations, the 
first Fourier component clearly dominates the response 
up to amplitudes close to the mean rate, where nonlin- 
earities are naturally expected, as there are no negative 
firing rates. This demonstrates that the linear response 
theory, although rigorously valid for small modulation 
amplitudes only, predicts the response in a large dynam- 
ical range of input amplitudes. 



IV. SUMMARY AND DISCUSSION 

The dynamical response properties of the generalized 
6'-neuron with adjustable AP onset speed were calculated 
in the presence of a fiuctuating correlated background 
noise. Methodologically we introduced a new approach 
which is based on the expansion of the corresponding 
Fokker-Planck operator to a complete set of orthonormal 
functions, leading to a sparse matrix representation. We 
then computed the eigenvalues and eigenfunctions of this 
matrix using an iterative scheme, the Arnoldi method. 
The high frequency limit was calculated analytically. It 
turned out, that the response amplitude decays as uj~'^ , 
where 7 depends on the kind of stimulation and, surpris- 
ingly, the phase at which a spike is emitted. As soon 
as this point differs from tt, where the dynamics is in- 
sensitive to external inputs, the exponent 7 is 1, giving 
the same asymptotic response behavior as conductance 
based neuron models. Using the eigenvalues and eigen- 
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Figure 4: Response amplitude for increasing values of the action potential onset speed. In the left column the response of the 
^-neuron for modulations in the mean input current and the noise amplitude is shown for different values of the stationary 
firing rate. The response exhibits a cut-off behavior, frequencies larger than the stationary firing rate are strongly damped. 
The middle column shows the response of the generalized model for both types of modulation and a stationary rate of 20Hz. 
For increasing values of the action potential onset speed the response amplitude grows for frequencies in the interval from 
lOOHz to IkHz, while the resonance maximum only slightly shifts to larger frequencies. The response of the noise modulation 
is much larger in this interval than the response to modulations in the mean current. For comparison the right column shows 
the response of the leaky integrate-and-fire (LIF) model. The response amplitude does not decay for large frequencies, for 
modulations of the noise amplitude it can even grow with increasing input frequencies. Parameters in the LIF simulation are 
as in _7j, except Ts = 10ms, a = 5mV, lo = 14.6; 16.2; 17.5; 19.5mV for a mean firing rate of 2, 5, 10, 20Hz. 



functions we then presented a method to evaluate the 
dynamic response to small time-varying inputs. There 
we found that for the classical 0-neuron model the re- 
sponse exhibits a cut-off behavior: For a modulation in 
the mean input current as well as for a modulation in the 
noise amplitude frequencies above the stationary rate of 
the neuron were strongly damped. In the generalized 
0-neuron the damping in the regime up to IkHz is sub- 
stantially reduced for both types of input modulations 
when the AP onset speed is increased, although the high 
frequency limit is the same as in the classical 0-neuron. 
The response amplitude for the noise amplitude modula- 
tion is typically much larger than the response amplitude 
for the mean input modulation. The Hnear response the- 
ory, although only derived for small modulations of the 
input current turned out to be valid in a large dynamical 
range, which we demonstrated by direct numerical sim- 
ulations. Amplitudes of the rate modulation up to the 
mean output rate turned out to be well described by the 
linear theory. 



Simple phenomenological, yet dynamically realistic 
models of cortical neurons are of key importance for 
studies in theoretical neuroscience, starting from stud- 
ies on spike timing to large scale network simulations or 
analytical network studies. While stationary response 
properties, such as mean firing rates or processes have 
been studied in many models, which operate on long 
time scales, e.g. adaptation (see e.g. il8, 33j), stud- 
ies on the dynamic response properties are rare. Most 
of these studies consider the dynamic response in the 
class of integrate-and-fire (IF) models 13, 21, 2^. In 
these studies, it was demonstrated that IF models can 
relay incoming stimuli instantaneously. Recently it was 
shown, however, that this response behavior strongly dis- 
agrees with the response of conductance based models 
and rather represents an oversimplification of the model 
than a feature of real neurons [l4.l2Sl|. While in [2^ the 
response properties of the classical 6'-neuron were inves- 
tigated, the authors of 0| studied another phenomeno- 
logical neuron model, the EIF model, which mimics the 
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Figure 5: Amplitude of the first four Fourier components as 
a function of overall modulation amplitude of the population 
averaged firing rate for (a) modulations in the mean input cur- 
rent and (b) for modulations in the noise amplitude {a = 0.7). 
The mean output rate is 20Hz, the modulation frequency IHz. 
The dashed line is the diagonal. Up to amplitudes close to 
the mean output rate, the first Fourier component is indistin- 
guishable from the diagonal, indicating that the response is 
essentially linear. Starting from amplitudes comparable to the 
mean rate, the infiuence of higher order Fourier components 
becomes substantial. The insets show the rate modulation for 
a modulation amplitude of 20Hz. 



dynamical response properties of a conductance based 
model. Our study corroborates and extends some of their 
results using a generalized model of the classical 6'-neuron 
0, E^l, a canonical model of conductance based neu- 
ron models, which exhibit type-I excitability and which, 
in contrast to IF models, incorporates a dynamic ac- 
tion potential onset. While the classical 6'-neuron model 
was ori gina lly studied in the super-threshold, noise-free 
case[i3,[l3|; recent studies focused on the resp onse in the 
presence of fluctuating input currents [allallal- These 



studies indicated that in a large parameter regime the 9- 
neuron exhibits the same stationary response properties 
as cortical neurons, e.g. a realistic shape of the f-I curve 
and irregular flring in the subthreshold regime. 

Despite these results, a major point of criticism ques- 
tioning the biological relevance of the model, remained: 
While the 0-neuron reflects the dynamics at the onset to 
repetitive flring, it lacks the sharp action potential up- 
stroke found in more detailed conductance based models 
and real neurons f3|- It was further argued that this 
deflciency results in a high frequency limit of the lin- 
ear response amplitude, which decays too fast cx 
while the linear response amplitude in conductance neu- 
ron models only decays oc uj~^. To address these issues 
we generalized the classical 6'-neuron, incorporating an 
adjustable AP onset speed, thereby mimicking the fast 
sodium activation at the action potential onset. Surpris- 
ingly, our study reveals that the high frequency limit, 
does not depend at all on the speed at the AP onset, but 
rather on the phase variable, at which action potentials 
are emitted. If at this point the dynamics is insensitive 
to external inputs, as in the classical 0-neuron, the de- 
cay of the linear response amplitude is at least oc 
whereas the decay is always oc uj~^ if the dynamics is 
not completely insensitive to external inputs, as is the 
case in conductance based neuron models. Moreover, the 
full transmission function reveals that the onset of the 
high-frequency limit can be shifted to very high frequen- 
cies if the speed of the AP onset is increased. These 
results question the relevance of the high frequency limit 
as a criterion for the typical transmission speed of neuron 
models. 

For the computation of the Hnear response amplitude 
we did not resort to direct numerical simulations, but 
used a method based on the eigenfunctions and eigenval- 
ues of the Fokker-Planck operator, describing the dynam- 
ics of the probability density function in the presence of 
a temporally correlated fluctuating input current. While 
this approach is in general well-known (see e.g. and 
[2^ for an appHcation to the non-leaky integrate- 
and-flre model in the presence of an uncorrelated back- 
ground noise), we derived a sparse matrix representation, 
for which we computed eigenvalues and eigenfunctions 
with very high numerical accuracy using a fast iterative 
scheme, the Arnoldi method QHS]- Compared to pre- 
vious studies on dynamical responses |7t ll^ll4j . this al- 
lowed for the computation of the linear response proper- 
ties with an accuracy that would be hard to meet by a 
direct simulation of the single neuron dynamics. 

Besides this, our results provide a direct link to exper- 
iments. In a recent study [j| it was shown that the AP 
width in neocortical neurons is strongly correlated with 
the critical frequency up to which a neuron can phase 
lock to sinusoidal input stimulations. This is indeed the 
same result we found for the generalized 6'-neuron: With 
increasing AP onset speed the response amplitude shifts 
to larger frequencies, enabhng the model to respond to 
frequencies much larger than its own stationary rate. In a 
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second experimental study it was demonstrated that cor- 
tical neurons subject to fluctuating input currents adapt 
their instantaneous firing much faster when stimulated 
with a step input in the noise amplitude than with a step 
mean input current j^jj- This behavior is well reproduced 
by the generaHzed 6'-neuron. For increasing values of the 
AP onset speed, the response amplitude at high frequen- 
cies is one order of magnitude larger for the stimulation 
in the noise amplitude, compared to the stimulation in 



the mean input current. Both results strongly suggest 
that the generaHzed 0-neuron, despite its simpHcity and 
analytic tractability, captures well the essence of the AP 
generating mechanism of multidimensional conductance 
based neuron models. Future experimental studies will 
have to show to what extent the generalized 6'-neuron 
predicts the dependence of the dynamical response prop- 
erties on the AP generating mechanism. 
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